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We develop an extension of the well-known BCS-theory 
to systems with trapped fermions. The theory fully includes 
the quantized energy levels in the trap. The key ingredient 
is to model the attractive interaction between two atoms by 
a pseudo-potential which leads to a well defined scattering 
problem and consequently a BCS-theory free of divergences. 
We present numerical results for the BCS critical temperature 
and the temperature dependence of the gap. They are used 
as a test of existing semi-classical approximations. 

Considerable interest in the field of ultracold gases has 
been sparked by the achievement of Bose-Einstein con- 
densation in the bosonic systems 87 Rb, 23 Na and 7 Li in 
1995 Q. Recently, several experimental groups have ex- 
tended these experiments to the case of trapped fermions. 
As a first step, it is attempted to achieve a degener- 
ate Fermi gas. In a possible next step, the celebrated 
Bardeen-Cooper-Schrieffer (BCS) phase transition could 
be observed. A promising candidate for achieving this 
transition is the isotope 6 Li: By trapping Li in two hy- 
pcrfine states, one can take advantage of the strong (at- 
tractive) interactions due to s— wave scattering between 
atoms in different hyperfine states. 

BCS-pairing occurs in a multitude of physical systems 
(e.g. electrons in metals, electron- hole exciton systems 
and neutron-proton systems) which share the character- 
istic that the formation of bound states (Cooper pairs) 
between the strongly coupled constituents is energetically 
favorable ||. For ultracold atoms in traps, the interac- 
tions are much better-known than in most of the above 
mentioned systems j3|. Therefore, the achievement of a 
superfluid state in these systems opens up the possibil- 
ity of testing our theoretical models, in particular, the 
validity of the BCS theory itself. Furthermore, the in- 
teraction strength and the density of the gas are exper- 
imentally tunable which, in principle, makes it possible 
to study the crossover from BCS pairing to Bose-Einstein 
condensation of bosonic pairs Q . 

To describe experiments with trapped fermions, one 
has to extend present theories (i) to include the discrete 
nature of the quantum energy levels of the particles in the 
trap, and (ii) to take into account the interactions spe- 
cific to the atomic case. We consider a model of trapped 
fermions with two internal states. At low energies, the 
p-wave interaction between atoms in the same internal 
state is negligible compared to the s-wave interaction be- 
tween atoms in different internal states. The latter in- 
teraction is characterized at low energies and for dilute 
gases (fci?r e < 1, where r e is the effective range of the 



interaction, kp the Fermi wavevector) by a single param- 
eter, the scattering length a. In order to achieve pair 
formation the interaction has to be strong; here we as- 
sume | a | 3> r e §. In this case, an excellent model for the 
atomic interactions is provided by the pseudo-potential 
discussed in ||. This model potential allows us to ob- 
tain an extension of the BCS theory to inhomogeneous 
systems which is free of divergences. 

We use this theory to calculate various observables in- 
cluding the critical BCS temperature and the density dis- 
tribution of the gas. We predict that this transition oc- 
curs at experimentally accessible densities and temper- 
atures, due to the large negative scattering length for 
6 Li atoms. We also carefully compare the results of this 
general theory to those of a theory based on the Thomas- 
Fermi approximation |?J (see also [||). In that way, we 
establish a region of validity of this approximation. 

Our model Hamiltonian for the trapped atomic gas 
includes only interactions between atoms in different in- 
ternal states: 

+ Jd 3 Rd 3 rd 3 r' (r'\V RM \r) 

x ^(R + ^U(R - ~)V>-«r(R - ^)W(R + \). (i) 

Here the fermion field operator Vf( r ) annihilates a 
fermion in the position eigenstate |r) with an inter- 
nal state a = +,— . The single particle Hamiltonian 
Ho = — l^-V 2 + Uo(r) — [i includes the trapping potential 
C/o(r). We assume there is an equal number of particles 
N a in each state and hence a single chemical potential 
/i |J. In Eq.(Q), we keep a general non-local interaction 
Vrm assuming only that it does not affect the center of 
mass motion R of the interacting particles. 

The interaction between atoms is often approximated 
in the center of mass frame by a contact potential, that 
is Vrm = 4:Ttah 2 S(i)/m. However, this approximation 
leads to an ultra-violet divergent theory. This reflects 
the fact that the contact interaction is an effective low- 
energy interaction invalid for high energies. One way to 
remove this divergence is to introduce an energy cut-off in 
the interaction. This approach has been used to describe 
superconductivity in metals where a natural cut-off in the 
form of the Debye frequency exists. Another method is 
to express the coupling constant in terms of the two-body 



1 



scattering matrix obtained from the Lippman-Schwinger 
equation. This renormalization scheme has been imple- 
mented in the literature only in the homogeneous case 
j|. Here, we put forward a technique valid also in the 
inhomogeneous case. We use the pseudo-potential Vrm 
p] defined for an arbitrary function (j>(r) by 



<r|Vfl A/ |0) = g6(r)d r [r<f>(r)} 



(2) 



with g = Airah 2 jm. We first note that in contrast to the 
usual contact potential the pseudo-potential leads to a 
well defined two-body scattering problem; the scattering 
state in the center of mass frame for two particles with 
a relative momentum p = (pi — P2)/2 and a relative 
position r is 



(r) 



Jp-r/h 



a ipr/h 



1 + ipa/h 
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For a potential of finite range r e the above form, diverging 
as 1 jr for r — > 0, is only valid for r»r e and for pr e /h <C 
1 The pseudo-potential has an effective range r e — 0; 
this does not lead to any physical problem as the 1/r 
divergence is regularized by the operator d r [r-] in Eq.(||). 

We follow the steady-state mean-field approach of BCS 
theory the quartic terms of H are replaced by 

quadratic terms taking into account all possible binary 
contractions in the spirit of Wick's theorem 0]. This 
leads to linear equations for the field operators, 

d f 

± J d 3 r'A(r,r')^(r'). (4) 
The Hartree field W(r, r') is defined as: 



W(r,r')= / d 3 y(r-r'+^|V iw |r'-r+^> 



x (4( r '-|)V; ± (r-|)) 



(5) 



(W is the same for both internal states). The pairing 
field describes correlations due to Cooper pairing: 



A(r,r') = / d 3 y(r-r'\V RM \y) 

x (tji 



Eq.(^,^|J^) form a non- linear self-consistent problem. We 
now use the pseudo-potential of Eq.(^); as it introduces 
a S(r), we have A(r, r') = S(r — r')A(R) and W(r, r') = 
5(r - r')W(R), with R = (r + r')/2, so that Eq.(g) be- 
comes local: 

ihj^ ± (K) = [Ho + W(R)]i>±(R) ± A(R)4(R). (7) 



The regularizing operator d r [r-\ on the right-hand side of 
Eq.(||) plays no role for the self-consistent Hartree field, 
which is given by W(R) = <7(-0t (R)Vv(R))- Its perti- 
nence becomes clear for the pairing field: 

A(R) ee - 9 lima r [r{^(R + ^_(R-^]. (8) 

Here the operator d r [r-] is necessary as the expectation 
value (ip+ip-) in Eq.(|J) diverges as 1/r for r — ► 0. To 
see this, we calculate from Eq.(Q) the time derivative of 
(tp+ij)-) which vanishes as the system is in a steady state: 
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+ A(R + ir)^l(R+|)^_(R-|)> 

+ A(R - ir)(^(R - ^+(R + |)) - A(R)*(r) (9) 

where U = Uq + W. The presence of 6(r) imposes a 1/r 
divergence in the pairing field (V 2 (l/r) = — 4w5(r)): 

V r 777 

(</>+(R + -ty_(R - -)> = — ^A(R) + F reg (H) + 0(r). 
I l Ann r 

(10) 

This 1/r behavior in the pairing field could actually be 
expected from the 1/r behavior of the two-body scatter- 
ing wavefunction Eq.(|^). Since the pseudo-potential re- 
moves this divergence, Eq.(||) yields A(R) = —gF reg (R.). 

We now compare the prediction of the present theory 
with the standard theory for a homogeneous system [0. 
In this case Uq ee 0, and W, A do not depend on the 
position R. The pairing field for a temperature T can 
then be expressed as an integral 
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/t) 2 ] 1 / 2 where jl =_Jt + W. This integral diverges for 



To calculate F reg we add and 



with f{E) = [cxp^/fcs^ + l]- 1 and£; k = [A 2 + (%^ 

r = as in Eq.©, 
subtract from (t^ + ?/>_) the integral of a function hav- 
ing the same large k behavior as the integrand, that is 
iAe ik '7[^-/t+ie], e -> 0+. The integral of this func- 
tion is G M (r)A/2 where G^(r) is the single free particle 
Green's function. The contribution of G^A/2 to F reg is 
easy to calculate as is known explicitly. The remain- 
ing integral (ip+ip-) — G M A/2 converges for r — > 0. Using 
l/(X + ie) = V{l/X) - m8{X), we finally obtain: 
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The equation A = —gF reg then coincides with the gap 
equation obtained by renormalizing via the Lippman- 
Schwinger equation Q. 

We now turn to the inhomogeneous case, for which we 
have to use a numerical approach. Here, we assume an 
isotropic harmonic trap Uo(r) = ^muj 2 r 2 . Following the 
Bogoliubov technique [0 , we expand the field operator in 
eigenmodes defined by the mode functions (u^^Vri) solv- 
ing the eigenvalue problem: 

£>„(R) = [H + W(R)K(R) + A(R)^(R) 
E n v v (R) = -[H Q + W(R)}v v (R) + A(R)it„(R). (13) 

Here, is an elementary excitation energy and rj stands 
for (n,l,m), i.e. the principal quantum number n, the 
angular momentum I and the azimuthal quantum number 
m. For thermal equilibrium, the pairing function is 



GMR+-)tf_(R--)> = 
( ^(R+£)^(R-I)[l-2/(^)]. 



(14) 



and a similar expression holds for W . Inspired by the 
homogeneous case, we extract from Eq.(|l4|) the term 
A(R)/2 times the single particle Green's function 



G M (R,r) = <R 



1 



R 



2 1 Ho ' 2 
G7 9 (R) - 



0(r). 



(15) 



The key point is that the 1/r diverging term can- 
cels in the difference between the pairing field and 
A(R)G Al (R, r)/2 so that we can take the limit r — > 0. 
Expressing the Green's function in the harmonic oscilla- 
tor basis {| </>?)}, we obtain 

A(R) = {^(R)^(R)[1 - 2f(E n )] 



A(R) |0°( R )I 2 
2 (n + §)?kj - n 



A(R)G^ (j R) (16) 



where we recall that 7/ = (n, I, m). In a practical numer- 
ical calculation, the infinite sum in Eq.(|l6|) is, of course, 
replaced by a finite one. The regular part of the Green's 
function can be expressed in terms of a one-dimensional 
integral which can be calculated numerically. 

Eq.(p"3|) for the (u^,u^)'s, the self-consistent deter- 
mination of A(R) given by Eq. ( |l6| ) and the one for 
W(R) constitute a non-linear problem. To obtain the 
critical temperature T c , one linearizes this problem for 
small A(R) which leads to the integral equation A(s) = 
J d 3 r M(s, r) A(r). The temperature for which the high- 
est eigenvalue of the kernel M crosses 1, is then T c [|. 

We have performed a numerical diagonalization of the 
kernel M for a varying chemical potential /i. For the 



interaction, we took the parameters of 6 Li that is a — 
— 1140 A jll| and a trapping frequency of 820 Hz which 
gives g — —l 3 hu> with I = (h/mui) 1 / 2 . In Fig|l], we plot 
T c as a function of N a . As hu/k B — 40 nK the calculated 
critical temperature seems experimentally obtainable. 

To obtain the spatial structure of the pairing field 
A(R) for arbitrary temperatures, we solve numerically 
the whole self-consistent non-linear problem. We plot 
m Fig|| the pairing field as function of R for a rela- 
tively large number of particles and a temperature much 
smaller than T c . In this low temperature regime the 
Cooper pairing takes place over the whole trapped cloud. 

In both of the above figures, we also compare with the 
Thomas- Fermi approximation (TEA), in which the sys- 
tem is treated as being locally homogeneous [0 , neglect- 
ing the discrete nature of the energy spectrum. There 
are two conditions for the validity of the TFA. First, 
the correlation length between the unpaired fermions 
~ 1/kp should be much shorter than the spatial radius 
ttf = y / 2fi/muj 2 of the cloud, requiring /i 3> Tllo. Also, 
the size £ of the Cooper pairs must be much smaller than 
ttf making a local theory for the pairing reasonable. 
For T = 0, kpt; ~ /i/A^-o ~ fi/ksT c where we have 
used A T=0 = 1.76k B T c § valid in the TFA. For T = T c , 
we have the same estimate kpt; ~ fj,/kBT c . For any T 
below T c the inequality £ <^.ttf then reduces to the re- 
quirement that the critical temperature must be much 
larger than the trap level spacing for the TFA to work; 
i.e. kBT c /Ttio ^> 1. 

We see from Fig.l and Fig. 2 that the agreement with 
the TFA is reasonably good for ksT c > fiui. To deter- 
mine the region of validity of the TFA more clearly, we 
plot in Fig. 3 the quantity S = J d 3 R A(R) as a function 
of temperature. In Fig|| (a), we chose a small number 
of trapped particles; in Fig.y (b) a much larger num- 
ber of particles are trapped. The temperature where the 
gap, that is S, vanishes, determines T c . For (a), we find 
k B T c = OAShuj, for (b) k B T c = 2.8huj. In Fig.3 (b), the 
agreement with the TFA is reasonably good. In Fig.3 
(a), the prediction of the TFA is larger by two orders of 
magnitude; this is not surprising given that the condition 
k B T c 3> Tilo is not satisfied. 

In conclusion, we have implemented a BCS theory for 
a dilute gas of weakly interacting fermionic atoms in a 
trap. It reproduces the well-known regularized gap equa- 
tion for the homogeneous case, and it provides a method 
for achieving a finite theory for a trapped gas taking into 
account the discrete nature of the normal state trap lev- 
els. Based on this theory, we intend in a next step to 
include time dependence in our treatment to study the 
response of the system to external perturbations, in a 
search for observable signatures of the BCS transition. 
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FIG. 2. Gap function A(R)//iw as function of R/l for 
g = -l 3 hu. We have fi = 31. 5hu (yielding N a ~ 8000) and 
fcgT = ksTc/10 = 0.28?kj. Solid line: numerical solution of 
the complete BCS theory. Dashed line: TFA. 
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FIG. 3. Size of the gap parameter in units of Tuol 3 as 
function of ksT/Tiuj for a scattering length g = —l 3 hui, and 
for (a) /i = WJahtd yielding N a — 300, and (b) u, = 31.5hu> 
yielding N a ~ 8000. Solid line: numerical solution of the 
complete BCS theory. Dotted line in (b): TFA. 



FIG. 1. The critical temperature kBT c /hw as a function 
of the number of particles N a in each of the internal states for 
g — —l 3 Ttio. The solid line is obtained from numerical solution 
of the linearized gap equation. The dashed line depicts the 
result of the TFA in the parameter space where it is valid (i.e. 
k B T c > huj). 
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Figure 3(a) 
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